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Abstract. 

We derive a trace formula that expresses the level density of chaotic many-body 
systems as a smooth term plus a sum over contributions associated to solutions of 
the nonlinear Schrodinger (or Gross-Pitaevski) equation. Our formula applies to 
bosonic systems with discretised positions, such as the Bose-Hubbard model, in the 
semiclassical limit as well as in the limit where the number of particles is taken to 
infinity. We use the trace formula to investigate the spectral statistics of these systems, 
by studying interference between solutions of the nonlinear Schrodinger equation. We 
show that in the limits taken the statistics of fully chaotic many-particle systems 
becomes universal and agrees with predictions from the Wigner-Dyson ensembles of 
random matrix theory. The conditions for Wigner-Dyson statistics involve a gap in the 
spectrum of the Frobenius-Perron operator, leaving the possibility of different statistics 
for systems with weaker chaotic properties. 
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Feynman’s path integral (see e.g. |T]) provides a convenient way to represent 
the propagator of a quantum mechanical system, and an excellent starting point 
for semiclassical and related approximations. Prime examples are van Vleck’s 
approximation of the propagator of a quantum system as a sum over contributions of 
classical trajectories [2], and Gutzwillcr’s seminal work [3] relating the energy spectrum 
of chaotic single-particle systems to periodic classical trajectories. These semiclassical 
methods provide one of the foundations of the field of quantum chaos ® 0 E]. For 
a many-particle system identifying the semiclassical limit is less obvious. A promising 
approach is to consider the path integral in second quantisation, running over different 
choices for the macroscopic wave function parameterised by position and time. One can 
show that in the semiclassical limit and in the limit where the number of particles N 
is taken to infinity this path integral is dominated by stationary points of the action, 
corresponding to solutions of the nonlinear Schrodinger equation or Gross-Pitaevski 
equation. These solutions take on a role analogous to the one played by classical 
trajectories in the single-particle theory. However previous work usually focused either 
on studying the full problem in second quantisation or, as e.g. in nuclear dynamics |7], 
on one dominating solution to the Gross-Pitaevski equation m- This does not exhaust 
the power of the approximation. In particular keeping the sum over different solutions 
of the Gross-Pitaevski equation allows to account for crucial interference effects between 
such solutions. An approach where the semiclassical propagator of bosonic many-particle 
systems was used to study these effects was pioneered in [10] for coherent backscattering, 
see also HQ for applications to fermionic systems. 

In the present paper we will focus on a further fundamental problem for which 
the interference between stationary points of second quantised path integrals is of vital 
importance: the statistics of the energy levels of many-body systems. To do so we will 
first derive an approximation of the level density in terms of stationary points of the 
action, and then study the interference between these points. To our knowledge the 
consequences of interference effects for many-body spectral statistics have not yet been 
investigated explicitly. We will see that this statistics depends crucially on the dynamics 
generated by the Gross-Pitaevski equation. If the dynamics is fully chaotic (in the sense 
to be specified below) the statistics in the limits considered becomes universal, and 
agrees with predictions from random-matrix theory (RMT). These predictions entail, 
for instance, a repulsion between the energy levels. 

This universal behaviour mirrors the behaviour of chaotic single-particle systems 
studied in the semiclassical limit. For such systems spectral statistics faithful to random 
matrix theory was conjectured in [12], and a semiclassical explanation was developed 
in da ei m na na [18]. This explanation is based on Gutzwiller’s trace formula 

[a 0 0 eo 

d(E) « d(E) + —Re V A p e iS ^ h 
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which expresses the level density as a smooth term d(E) plus fluctuations associated 
to classical periodic orbits p (with energy E ) of the system. Here S p is the reduced 
action of the orbit given by S p = j p(t) ■ q(t)dt where q[t) denotes the vector of 
generalised coordinates and pit ) denotes the associated momentum. The amplitude 

yprim -iMpf 

A p = —^===^=j- depends on the primitive period TP nm , the so-called Maslov index p p 

and the stability matrix M p relating deviations in the end of the orbit p to deviations 
in the beginning; I is a unit matrix. (Throughout this paper I will denote unit matrices 
with a subscript indicating their size if it is not clear from context.) Our first aim will be 
to generalise the trace formula to bosonic many-particle systems in second quantisation, 
with solutions of the Gross-Pitaevski equation taking the role of classical trajectories. 

We will then use this result to investigate spectral statistics. An observable 
characterising spectral statistics is the two-point correlation function of the level density 
d(E) = S(E — Ej ) (where Ej are the energy levels of the system). This correlation 
function is defined by 



where the average (...) is taken over an interval of E for which d can be taken constant 
as well as over a small range of e. Inserting the trace formula one obtains a double sum 
over solutions of the Gross-Pitaevski equation, and by taking into account interference 
between solutions we indeed recover statistics in agreement with RMT. 

More precisely this agreement holds for the statistics inside appropriate subspectra 
defined by the symmetries of the problem; for many-body systems we have at least one 
symmetry, particle number conservation, requiring to consider subspectra associated to 
a fixed particle number. Further refinements (to be discussed below) arise in case of 
geometrical symmetries. The precise ensemble to be chosen depends on the behaviour 
of the system under time reversal. The most frequent case involves systems invariant 
under a time-reversal operator squaring to 1. In this case (assuming that there are no 
further symmetries) one has to use Wigner’s and Dyson’s Gaussian Orthogonal Ensemble 
(GOE), i.e. predictions for spectral statistics are obtained by modelling the Hamiltonian 
through a real symmetric matrix, and then averaging over all possible such matrices 
with a Gaussian weight. In the absence of time-reversal invariance one averages instead 
over the ensemble of Hcrmitian matrices with a Gaussian weight, the Gaussian Unitary 
Ensemble (GUE). 

A paradigmatic example for the systems to be considered is the Bose-Hubbard 
model, a model with L discrete sites (labelled by k — 0... L — 1) accommodating bosonic 
particles. Denoting the creation and annihilation operators for particles at these sites 
by a\ and a± the second-quantised Hamiltonian can be written as 

H = ~^ y^(Qfc + jQfc + alhfc+i) + y ) 2 d 2 k (3) 

k k 

describing hopping between the sites as well as interaction between particles on the same 
site. One can consider periodic boundary conditions and set cil = a 0 , a) L = aj. More 
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generally we are interested in discrete bosonic many-body Hamiltonians that are of the 
form 

H ^ ^ hkici k di T ^ ] Uki mn d k did m d n (4) 

kl klmn 

or have even higher-order interactions; here real coefficients imply time-reversal 
invariance. 

A crucial requirement is that the underlying ‘classical dynamics’ is chaotic. This 
dynamics is obtained by replacing the creation and annihilation operators by mutually 
complex conjugate time dependent variables ifl, ifk where if = (V’o, ifi, • • •) can be 
interpreted as a macroscopic wave function and ifk as its value at site k. As N = ^2 k a\dk 
is the particle number operator, the macroscopic wave function is normalised to have 
Sfc \^kY = N where N is the particle number. The associated analogue of Hamilton’s 
equations for if k can be shown to read ihifk = jnY and takes the role of a discrete 
nonlinear Schrodinger equation. 

We note that the ‘classical’ Hamiltonian entering into this equation depends on the 
ordering of operators. For the Bose-Hubbard model the normal-ordered Hamiltonian in 
(J 3 J) yields an interaction term j |f/k| 4 ; however if the Hamiltonian is first brought to 
Weyl-ordered form (with all possible orderings of operators in a product contributing 
symmetrically) before replacing the operators by macroscopic wave functions one obtains 

fE,(l^l 4 -2|^l 2 + |)- 

For the Bose-Hubbard model the dynamics generated by the discrete nonlinear 
Schrodinger equation has been found to be mainly chaotic (with chaotic regions of phase 
space dominating compared to regular ones) in the case of several sites and comparable 
hopping and interaction terms [20]. In the same regime, numerical studies suggest 
spectral statistics in line with the GOE [2Tj, 22]; see also [23] for fermionic systems. 

To explain the observed faithfulness to RMT, we follow a semiclassical approach 
inspired by single-body spectral statistics as well as HD). Our first aim is to derive 
a trace formula for second quantised Bose-Hubbard-like systems. This will be done 
in Sect. 2, after a brief reminder of the corresponding derivation for one-body chaotic 
systems. Special emphasis is placed on the treatment of conserved quantities. In Sect. 3 
it is demonstrated how the obtained trace formula can be used in order to predict the 
spectral statistics of the system, especially the 2—point correlation function. We explain 
in detail how to generalise the approach previously used for one-body chaotic systems 
in Erst quantisation, and we discuss the range of validity of this approach. In Sect. 4 the 
consequences of discrete symmetries of chaotic many-body system, in particular for Bose 
Hubbard model are investigated. In Sect. 5 we present numerical results supporting our 
claims. Some more technical details of the derivation are put in two appendices. 


Spectral statistics of chaotic many-body systems 

2. Trace formula 


5 


2.1. Trace formula for single-body systems 

In order to prepare our derivation of a trace formula for bosonic many-particle systems, 
we want to briefly review the calculation leading to the trace formula for single-body 
systems. We refer to H U O EH EE] for further details. 

The level density can be accessed from the trace of the time evolution operator 

_ i_ fjf 

e fi via 

- i r°° 

d(E) = tr 5{E — H) = —Imi / dt e^ E+l0 ^ tr e~^ Ht , (5) 

7r h Jo 

where an infinitesimal positive imaginary part has been added to the energy to ensure 
convergence. The trace of the time evolution operator can itself be expressed as a path 
integral over phase space trajectories ( q(t),p(t )) 

tr e~* 6t = J D[q,p\e iRlq ’ p]/h , (6) 


where the action weighting each path is determined by the classical Hamiltonian H as 

R[q,P\ = f dt'(p(t') -q(t') - H(q(t'),p(t'))). (7) 

Jo 

Now we are interested in an approximation of this path integral in the semiclassical 
limit, i.e., the limit of large quantum numbers implying that typical classical actions are 
much larger than H: formally this limit is often denoted by h —> 0. In the semiclasssical 
limit the path integral is dominated by stationary points of the action, corresponding 
to periodic orbits that satisfy Hamilton’s equations of motion, and a stationary-phase 
approximation leads to a discrete sum over such periodic orbits. For chaotic dynamics, 
free from continuous symmetries, the periodic orbits are isolated but one has to take into 
account that each of them formally gives rise to a one-parameter family of stationary 
points distinguished by which phase-space point along the orbit is taken as the initial 
one associated to t = 0. The Laplace transform in ([5]) can then be carried out in a 
further stationary-phase approximation and one obtains the trace formula 


d(E) 


« d(E) + 



yprim g-i/Xpf 


v'l det(AT p -1)| 


e iS p /h 


( 8 ) 


already shown earlier. Here the determinant and the phase factor involving the 
Maslov index arise from the Hessian matrix of the action entering the stationary-phase 
approximation, and the factor T£ rim results from integration over different choices of 
initial points. The summand d(E) can be derived from a careful treatment of the lower 
limit of the time integral (|5]). 
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We now want to generalise the trace formula to many-particle systems in second 
quantisation, for the case of bosonic particles at discrete sites. In second quantisation 
it is natural to work in a basis of coherent states, i.e., normalised joint eigenstates of all 
annihilators dk with eigenvalues fjk- Further differences from the single-particle setting 
arise from operator ordering as well as the conservation of the particle number. 

Again we access the level density from the trace of the time evolution operator 
using (5). For many-particle systems the latter trace is given by the path integral 




tr e *- H ' = / D[xjr, *}e 


iR[4>,g>*]/h 


(9) 


over all macroscopic wave functions ■*/?( t ') that return to their initial value after time t. 
Here the action is 

^ dt' 0 '^(0 - . (10) 

Eq. (J9]) is related to but somewhat simpler than the path integral for matrix elements of 
coherent state time evolution operator EH; our use of tre * Ht is motivated by [2l|25]. 

In case of normal ordering the path integral can be derived by splitting the time 
interval t into J steps of width r = f- and then using the result for the short-time 
propagator 

(4> j+ i\e~ lHT/h \il>j) = exp 


</b*+1 • ~ 


l^i+il 2 + W 2 i jj , ,* /N 

——2— “ n H ^ j+ i’^ )r 


+ 0(r 2 ); 


( 11 ) 

after integration over the macroscopic wave functions at all time steps this leads to a 
discrete path integral with the action 

R[^A*} = Y ( _ 7^*+i ■ (^i+1 _ ^i) _ (12) 


and then to (10) after taking the limit J —> oo. To £x notation, throughout the paper 


j will be a time index, k will label degrees of freedom e.g. associated to sites, and bold 
vectors assemble all choices for k. j is defined modulo J and k is taken modulo L. With 
these conventions the integration measure in the path integral is given by Y\j k • 

The integral thus obtained agrees with the phase-space formulation of the path 
integral in first quantisation if the macroscopic wave function is related to canonical 

coordinates and momenta through + Wk) or ^ = \Jdfe~ l9k . We can now 

perform a semiclassical approximation valid in the limit H —>■ 0. Here is taken of 


the order such that the first term in the action (10) becomes independent of h; 
the coefficients in the Hamiltonian are scaled following J/h t —> J, U/h 2 i —> U in such 
a way that the Hamiltonian is independent of h as well (We will later also give an 
interpretation in terms of the limit N —» oo.) In our limit the integral is dominated by 
periodic solutions of the nonlinear Schrodinger equation ihif k = (and its complex 
conjugate) following from the stationarity of R. 
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Importantly, however, our Hamiltonian has a continuous (gauge) symmetry w.r.t. 
multiplying all components of the macroscopic wave function with the same phase factor. 
A consequence of this symmetry is that as mentioned the total particle number N is a 
conserved quantity commuting with the Hamiltonian. It is thus preferable to consider 
the density of levels forming the subspectrum associated to a fixed particle number. 

To implement this restriction we subject the Ik, Ok defined above to a canonical 
transformation. (See [ 281 [9] for alternative approaches.) This transformation is chosen 
to lead to I' k ,0' k where I' Q = J2 k Ik = ^SfclVk;| 2 = the remaining I k are linear 
combinations of the Ij, and the 0' k are defined such that the overall transformation 
becomes canonical. If we choose the transformation in such a way that the range of 

possible 0' k is limited to 27 t it is also convenient to let 0(. = As I' 0 is a 

conserved quantity the corresponding canonical coordinate 0' o must be absent from the 
Hamiltonian. The remaining variables with k > 1 parameterise a reduced phase space 
associated to the particle number N. 

We are interested in the part of the spectrum associated to a given value N of the 
particle number. The associated level density can be formally written as 

dx{E) = tr Sjy N 5{E — H) = — Imi / dt e ft(- B +»°) t tr 5^ N e~*- Ht (13) 

J 0 


with the Kronecker delta 


r-2?r 


5 *’ n ~ 2tt 


dtp e 


i<t>(N-N) 


(14) 


As one might expect, d]y(E) is determined by solutions periodic in the reduced phase 
space. To see this formally we split the exponential e^ N from (14) into factors e ir * >N ^ 
to be inserted into each of the short-time propagators 0- This gives an additional 
term 

</>*] = y ( 15 ) 

3 

to be added to the action, leading to 

dt'^^it'f 2 

in the limit J —y oo. The phase now becomes stationary if 


(16) 


dH ho , 

= wr 

where the addition only affects the dynamics of 9' 0 , replacing it by 

q, = dH _ 0 


dl’ 


t 


(17) 


(18) 


Hence, as anticipated, the requirement of periodicity is nontrivial only for the variables 
ip' k with k > 1. Jq is conserved and hence trivially periodic, and 0' o is required to be 
periodic only after modifying the dynamics through the ^-dependent term in the action, 
but all possible choices of 0 are integrated over. 
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We now have to determine the weight associated to each periodic solution. We recall 
that if the stationary points p of a given action R[x] (with x G M n for now) are isolated 
and we are taking the limit h —> 0, the integral over e lR ^ h can be approximated by the 
following sum over contributions associated to stationary points, 


d n x e iR[x]/h » ^(2tH) 


n/2 


detU^ 


h dx 2 


-1/2 


iRp/h if-ip 2 . 


(19) 


d 2 Rr. 


here p p is the number of the negative eigenvalues of the Hessian matrix for the 
stationary point p. If the stationary points are not isolated this leads to vanishing 
eigenvalues of the Hessian. In this case the integral over the directions orthogonal 


to the stationary-point manifold can still be computed using (19), but it has to be 
accompanied by an integral over the manifold itself. 


The stationary points of (12) are not isolated. In particular continuous time shifts 


of a solution of the nonlinear Schrodinger equation lead to a different solution, and the 
same applies to simultaneous shifts of all 0'- o . This is a consequence of the two conserved 
quantities, for the energy conjugate to time and the particle number conjugate to the 
However as a first step it is still helpful to compute the (rescaled) Hessian involving 
derivatives w.r.t. all components of the macroscopic wave functions at all time steps. 
Adopting a complex notation we consider 


H = \ 


d 2 R 


ftdW’o.V’o.V’I.i/’i.'") 2 ' 


( 20 ) 


Using the discretised action (12) the derivatives are given by 

d 2 R_ 


d *!) 2 Oil) 2 

d 2 R d 2 H(i/)*,il) 1 _ i) 


diJ) J 

d 2 R 


*2 


dij) 


*2 


-r, 


ft, 9 2 »(«/>y 1,^-) 

= -I-T 


di/)* j+1 dil)j i di\)] +x di\) 3 

d 2 R __fr 

dij)*dil) j i 


( 21 ) 


where I = is a unit matrix of dimension L. The Hessian w.r.t. these variables then 
assumes periodic block tridiagonal form 

/ Ao Bq Cq \ 

c\ 

B n -2 

\B n —1 C n —i A n —\ J 


n = 


( 22 ) 
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A 2j - 


A 2 j+i - 




dip * 2 h 



+ 0(r 2 ) 
+ 0(r 2 ) 


B 2j = il 


B 2 j + \ — —il — 
C 2 j = —il — 


d 2 H(iP*,iP 3 ) T 

dip j dip* h 

dip* dip 3 h 


+ 0(t 2 ) 
+ 0(r 2 ) 


Cgj'+l — *1- 


(23) 


Here we neglected corrections of order r 2 arising from the fact that the arguments of 
the Hamiltonian in (21) are taken at slightly different times. We can now use a general 
formula for determinants of block tridiagonal matrices as in (22) that was derived in 
using a transfer matrix approach, 


det H = det(M - I 21 ) detip) 


(24) 


Here we have 


P — Bq ... B n _i 

?-i 


M = 


~B n _iA n _i — 5 n _ 1 C' n _i 
I 0 


-BZ'A„ -n-'c, 


(25) 


In Eqs. (24), (25) we have slightly modified the numbering of indices from [29] and 


taken the dimension of the block matrices as L. Due to n — 2 J the sign factor in (24) 


is just equal to 1. To evaluate M we group the factors in (25) into pairs 


M - I ~ B 2j+l A 2j+l ~B 2j 1 +1 C 2j+ i 

1 “ ‘ I 0 


—B 2 jA 2j —B 2 j 1 C 2j 


(26) 


which using (23) can be simplified to 


IT 

Mj = I 2 L + -jT 


dip j dip* 
d^iip*,^) 

dr* 


9 2 H(iP*,ip j ) 

dipp 

d 2 H(ip*,ip j ) 
dip j dip j 


+ 0( r 2 ). 


(27) 


One can now show that multiplication with Mj maps small deviations (dip*, dip j) from a 
given solution of the nonlinear Schrodinger equation at time jr to the resulting deviation 
at time (j + l)r. This follows immediately by linearizing around the equation 


ipj +1 = '•Pj + + °( T2 ) = tpj 


dH ir 
dip* h 


+ 0(r 2 ) 


(28) 
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M — ... M\Mq 


( 29 ) 


(understood in the limit J —y oo) maps deviations in at time 0 to those at time t. 
The matrix P is given by 

t-,-/ ir<9 2 P(0*,0 7 ) , 

P = Bq ... B n _\ = n 1 - IT + 0(r 2 ) 


= IJ exp 

3 




h dip'dipj 

+ 0(r 2 )V 


h d'4>*d'4> j 7 

evaluating its determinant and taking the continuum limit then leads to 

( 9 - 0 * ( 9-0 


detP = exp ( ——tr I dt 


(30) 


(31) 


The factor (detP) -1 / 2 arising from this term in (detP) -1 / 2 is known as the Solari- 
Kochetov (SK) phase; it is in line with previous work about the propagator in normal 
ordering [50, .5TI \T7\. 


2.5. Treatment of the conserved quantities 

We now modify our treatment to take into account particle number and energy 
conservation. 

For notational convenience it is helpful to complement our earlier canonical 
transformation singling out the particle number by a further transformation affecting 
only the variables with k > 1. This transformation is defined only in the vicinity 
of a periodic solution and leads to I[ indicating the energy on the orbit, and 6\ the 
time along the orbit. The remaining variables with k > 2 indicate transverse 

deviations from this orbit. If desired they can again be turned into complex variables 
0(. as abov^f and the vector formed by all 0[. with k > 2 will be denoted by 0^. The 
present transformation is analogous to the transformation to parallel and perpendicular 
coordinates in the derivation of the trace formula in first quantisation [31 31EJ- We note 
that it cannot be expanded to a transformation in the full phase space as this would 
imply integrability. 

When determining the weight associated to periodic solutions, it is now crucial 
to take into account that shifting all time coordinates 6' ]{ by the same amount leads 
to a valid solution, and the same applies to coordinated changes of the variables 6' ](] 
conjugate to the particle number. Hence there are two linearly independent ways in 
which continuous changes from one stationary point of the action lead to a different 

|| This requires suitable rescaling to make both quantities dimensionless and restrict the range of 0' k to 
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one. As second derivatives of the action in the associated directions must necessarily be 
zero the matrix % defined above must have a two-fold eigenvalue zero. 

To compare this to the behaviour of M — I 2 r we first of all observe that M maps 
deviations of the variables associated to k — 0,1 only to deviations associated to the 
same k. Written in terms of 9' k , I' k the stability matrix multiplies deviations ( 5Q' k , 5I k ) 
with 

1 b k 
0 1 


M<‘> = 


(32) 


i ; dAe k 

where b k = dI , ■ 


Here the diagonal elements indicate that the conserved quantities I k 
stay fixed and changes of 9' k are translated into equal changes in the end. In the right 
upper element A 9' k indicates the increase of 9' k along the orbit, e.g. the period for k — 1. 
The coefficient b k takes into account that a change of the energy typically changes the 
period of the orbit, and similarly a change of the particle_number typically changes the 

b k 


difference of the initial and final 9' 0 . As a consequence of (32), M^ — I 2 = 


0 


has 


determinant zero. 

To deal with these zeroes one can consider a perturbation of the Hamiltonian 
[25] that replaces one of the zero eigenvalues of the Hessian by a small value e. In 
the corresponding M^ a factor ej then enters in the lower left corner, turning the 
determinant into —ejb k . A brief account of these perturbative results for the present 
case is given in Appendix A The determinant of the Hessian matrix "H omitting 
the directions associated to conserved quantities can now be evaluated by considering 
perturbations for both k = 0 and k = 1 and dividing out the two factors e, leading to 


detTf = det(M — I 2 l- 4) det(P) l 7 2 6 0 ^i 


(33) 


Here M is defined in analogy to (27) and (29) but only w.r.t. the variables 
omitting k = 0,1. M maps initial deviations of 'ip ± *,'ij) ± to the corresponding final 
values. This meaning is precisely equivalent to the stability matrix appearing in the 
conventional trace formula. 

Our result for det 9i allows to evaluate (in a stationary phase approximation) the 
integral over all momenta and coordinates apart from 6*' 0 , 0' 1; as well as the fluctuations 
of 9b 0 , 0'j as j is varied. 

It remains to consider the constant (in j ) Fourier modes of Mi- Importantly, 
if we perform a discrete Fourier transform of 9’- k and want the associated Jacobian 
determinant to be 1, the integration variable parameterising the constant mode has to 

q/ 

be chosen as , i.e. \fj times the average of the 9’ jk . More natural parameterisations 
of the constant Fourier modes, such as by the average of the 0' k , entail a Jacobian \fj 
for each k = 0,1. These jointly cancel the from (detP)^ 1 / 2 . 

Integration over the constant modes now leads to multiplication with the integration 
ranges. For 9' 0 j the range is 27T, and for the time coordinates the range is normally the 
period. However if the periodic solution consists of several repetitions of a shorter 
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‘primitive’ periodic solutions, all distinctive stationary points of the action can be 
accessed by integration over the period T prim of the primitive solution only. Equating 
T p = T pnm if our solution does not consist of repetitions of a shorter one we thus obtain 
a factor T prim in all cases. 


We still have to evaluate the integral over 0 

— [ df e iR / h ~^ N . 
2n I 


(34) 


where the R from (16) equals 0/0 Due to = I' 0 the stationary phase condition gives 


Jo = hN as expected. Now due to (18) a change of 0 is equivalent to changing the 


obtained cancels with the one from (33). 


range A9' 0 by an opposite amount. This implies that = 


din 


dAd'g 


— —ft -1 

— 0 Q . 


The bn thus 


Altogether the trace of the time evolution operator, restricted to a given value of 
the particle number, is thus approximated by the following sum over solutions p of the 
nonlinear Schrodinger equation that are periodic with period T p = t in our reduced 
phase space 


H N e 




E 


'J'primpRp/H— 


(35) 


p v /27r ^l & ill det(Mp - I)| 

Here various factors from the integration measure, Eq. (19), Eq. (14), and the 6' 0 


integral have cancelled mutually. The result is in line with [.24>12S!- The Maslov index 
counts the negative eigenvalues of the part of the Hessian matrix R associated to 
k = 2,3,... (when brought to symmetric form involving derivatives w.r.t. /j fc , 9’ jk ). The 
index takes a similar role for k — 1; it is equal to 1 if b k < 0 and 0 otherwise. An 
analogous phase associated to k = 0 has already been cancelled by the stationary-phase 
approximation of (34). For notational convenience we also use the Maslov index to 


absorb the Solari-Kochetov phase jSDHSDEI! 


,7T 


/-4 o = -^ tr 


/o dil>*dil) 


(36) 


p 2 2 h 

We note that like the stability matrix also the action can be taken within our reduced 
phase space. The only difference between the two is the k = 0 contribution to the first 


term in (10), which is given by 

rt 


dt' ( —^ ip'* 0 (t')ip' 0 (t')] = [ dt'I' 0 (t')9' 0 (t') = 27 rhN- 


(37) 

>0 \ 0 /JO 

however for integer N this term has no influence on the phase factor e* R/ 4 Moreover, 
as the Hamiltonian is independent of 9j 0 , the change in the dynamics of 9' ]() mentioned 
above does not affect any of the quantities entering the trace formula. 


2.6. Level density 


Finally the Laplace transform of (35) 


d N (E) = tr 5pf N 8(E — H) — -^Imi J 


dt tv 5 N e * 


-iHt 


(38) 
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can be performed in a further stationary-phase approximation analogous to [3,0, .6] (and 
in fact similar to the integral above). The phase in e d Et +Rp)/ h becomes stationary if 


E = However as shown e.g. in [4] — is precisely the energy of the solution with 

period t. Hence the result of the stationary phase approximation will be a sum over all 
solutions that are periodic in the sense above and have fixed energy E rather than fixed 
period. For these solutions Et then cancels with the second term in the action © , hence 
the associated phase will be determined only by the first term S p = — | dt'xp*{t')-i^{t') 


dR 


also referred to as the reduced action. 

. (!) 


Given that (^-4 = 


-TT = —b x the factor 


\J‘lT\ih\b\ \ e l,J " p 2 arising from the stationary-phase approximation combines nicely with 


the one from (35). Altogether we thus obtain the anticipated trace formula 


1 


J^pnmj 


d]y(E) « d^iE) -f--ReV^ ._ 

A \/l det(A/ p — 


0 iS p /h 


(39) 


“v— 

— •An 


The summand dj^(E) gives the smooth part of the level density and arises from the lower 
limit 0 of the time integral. It is proportional to the volume of the energy shell in our 
reduced phase space and given by the pertinent variant of the Weyl or Thomas-Fermi 
formula 


ME) 


1 

(27 \h) L ~ l 


dl[... dI' L _ x de\... d0 L _ x 8{E - if (!', 6')). 


(40) 


This result (which is more meaningful if one avoids the canonical transformation in 
the beginning of subsection 2.5) is readily obtained using that for small times no time 
slicing is necessary. The action (12) then boils down to — H (t/?q, i/> 0 )t and the result 
follows directly after Laplace transformation. For the semiclassical approach to spectral 
statistics it will not be necessary to evaluate the - often involved - integral for diy(E). 
As in first quantisation [32] subleading corrections to this result are to be expected. 


2.1. Discussion 


Noting that \x^\ 2 = N the present approximation is valid not only for h —> 0 but also 
in the limit N —> oo. In the latter case it is helpful to rescale ?/? —>■ a fNij) to avoid 
changing the normalisation of the variables in the limit taken. If we want the hopping 
and interaction terms in a Bose-Hubbard like system to remain comparable we then have 
to adjust the corresponding coefficients in a way similar to the case h —> 0. This now 
entails J —> J, NU —> U. In case of the Bose-Hubbard model the resulting Hamiltonian 
entering the trace formula then satisfies 


= N 




(41) 


k k 

with U and J constant, the normalisation l^l 2 = 1 and the large parameter N and the 
whole action (12) being proportional to N. 
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If the Hamiltonian is written in Weyl instead of normal ordering one obtains an 
analogous result however without a Solari-Kochetov phase. This is in line with [24] as 
well as results for the propagator in [27J [33]. A derivation along the lines followed here 
will be given in Appendix B; it uses the analogue of the discretised action stated in 


and the corresponding Hessian is again evaluated with the help of 

An alternative derivation of the trace formula can be based on the semiclassical 
approximation of matrix elements (i/j^e lHt / h \^ W) in a basis of coherent 

states. Interestingly, in this case 'ipit') and i^*(t') first have to be treated as independent 
functions subject to the conditions i/>(0) = = only. The two 

functions become complex conjugate after evaluating the trace in a stationary-phase 
approximation, and the resulting trace formula coincides with the one obtained above. 


3. Spectral statistics 

We are now equipped to study spectral statistics. Inserting the trace formula into 
the definition of the two-point correlation function one obtains, as for single-particle 
systems, the double sum 

«(e) « 1 + Re ly , (42) 

w (,rM) / 

The correlation function is thus expressed in terms of periodic solutions of the nonlinear 
Schrodinger equation in a way that allows to keep track of crucial interference effects. 
The double sum can be evaluated in the same way as for chaotic single-particle systems. 

We now want to discuss this evaluation in more detail. In doing so we will emphasise 
the ingredients entering the calculation and check the conditions under which the 
reasoning for single-particle systems carries through to many-particle systems in second 
quantisation. For the details of the calculation for single-particle systems based on these 
ingredients we refer to the original literature quoted below as well as mm- 

3.1. Conditions 

The phase space of predominantly chaotic many-particle systems typically still has 
small stability islands. Hence it is important to stress that our theory describes the 
behaviour of states supported by the chaotic part of phase space. The spectral statistics 
is dominated by this contribution if the regular parts of phase space are small in 
comparison, as in the case of the Bose-Hubbard model in the regimes considered. 

In the chaotic part of the phase space our treatment requires a gap in the spectrum 
of the Frobenius-Perron operator. This condition implies various weaker requirements 
such as ergodicity and hyperbolicity. The Frobenius-Perron operator describes the time 
evolution of classical phase space densities [I], leading from a density po(x) at time 0 
to the density pt(x) = J d n x' 5(x — & t (x'))p 0 (x') at time t. Here ^(cc') gives the image 
of the phase space point x' under classical time evolution over time t. The leading 
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eigenvalue of this operator is 1, with the associated eigenfunction corresponding to a 
uniform density on the energy shell. The remaining eigenvalues can be written as e~ lmt 
where 7 m with Rey™ > 0 are the Ruclle-Pollicott resonances. The system is said to have 
a spectral gap if the remaining Re y m are bounded away from zero so that all modes 
associated to non-uniform phase space densities decay in time with at least a minimal 
rate. For many-particle systems the dynamics required to satisfy this condition is the 
one induced by the discrete nonlinear Schrodinger equation in the reduced phase space 
parameterised by I' k ,9' k (k — 1, 2,..., L — 1). 

As a further condition we assume for now that our system has no further symmetries 
beyond the particle number conservation already taken into account, however we will 
extend our treatment to deal with discrete geometric symmetries at a later point. 


3.2. Diagonal approximation 


When evaluating (42) it is important to look for pairs of solutions whose (reduced) 


actions S p , S p ’ are similar as this means that their contributions can interference 


constructively. In contrast, terms in (42) arising from pairs of orbits with large action 


differences oscillate rapidly as the energy is varied and are washed out by the energy 
average. 

The simplest pairs of solutions with similar actions involve two solutions that are 
identical (apart from the slight difference due to the offset in their energy arguments). 
If we neglect the difference of the corresponding factors A p , A p t and Taylor expand the 

dS 

exponent using that gives the period T p , the contribution from identical solutions 
can be written as the single sum 


-Rdiag(c) R® l ^ ' 


14 


[nhdy 


0 iT p e/ird 


(43) 


Under the condition of a spectral gap such sums over periodic orbits can be evaluated 
using a sum rule derived in the quantum chaos context by Hannay and Ozorio de 
Almeida In the notation used here it can be written as 




' o 


00 dT 
~T 


(44) 


where the dots represent an arbitrary property of the solutions that depends only on 
their period T. This rule is a general statistical property of periodic solutions in systems 
with a spectral gap; it is very helpful to extract information from these solutions even in 


situations where it is difficult to determine the solutions individually. Eq. (44) implies 


that even very long orbits give important collective contributions: while the factors \A p \ 2 
associated to these orbits decrease with increasing period their number increases, and 


both effects approximately compensate. Using (44) the sum in (43) can be evaluated to 
give 


Adiag (b ) 


1 

24 


(45) 
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Figure 1 . Pairs of orbits differing in encounters: Simplified sketch of a Sieber-Richter 
pair of orbits differing by their connections in a single encounter of two orbit parts, 
and a pair of orbits differing in two encounters involving two and three orbit parts. 
(Pictures from [17] (c )American Physical Society). 


This result, originally derived in [13], IT], is known as the diagonal approximation. It 
gives the first nontrivial term in an expansion of the two-point correlation function in 
/, after the leading term 1 present in (42). For time-reversal invariant systems we also 
have to keep track of mutually time reversed solutions, leading to a doubling of this 
result. 


3.3. Encounters 

Contributions of higher order in 1/e arise from pairs of orbits differing noticeably only in 
so-called encounters PI El. Inside these encounters two or more parts of the same orbit 
come close up to time reversal, and a partner orbit can then obtained by connecting 
the ends of these orbit parts differently and a subsequent adjustment to obtain a valid 
periodic orbit satisfying the equations of motion. For time-reversal invariant systems one 
also has to include the case where parts of an orbit are almost mutually time reversed. 
Two simplified sketches of such pairs of orbits are shown in Fig. [l] 

The existence of such pairs of orbits requires hyperbolicity which follows from the 
existence of a spectral gap. In a hyperbolic system the possible directions at each point 
in phase space (apart from the direction of the flow and the direction of increasing 
energy) are spanned by pairs of stable and unstable directions. Deviations between 
(parts of) trajectories along the stable directions decrease exponentially in time, whereas 
deviations along stable directions increase exponentially but decrease for large negative 
times. This allows to ‘change connections’ inside an encounter, by constructing a part 
of an orbit p' that moves away from one part of p and towards a different part of p\ its 
deviation from the first part has to be unstable whereas the deviation from the second 
part is stable. 

When deriving the contribution of such pairs of orbits for single-particle systems, 
the deviations between encountering parts of an orbit were measured in a system of 
coordinates associated to the stable and unstable directions [351 mi [37]. This carries 
over in the present scenario as well. For example, if two parts come close we consider the 
points where these two parts pierce through a Poincare surface of section in our reduced 
phase space. By transformation from the difference between the /(,, 6' k (k = 2,..., L — 1) 
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we then obtain L — 2 pairs of coordinates Sk,u k characterising the deviations between 
the two orbit parts in the stable and unstable directions. If the parts are almost time- 
reversed instead of close in phase space we instead have to consider the deviation of one 
part from the time-reversed of the other. 

We can now determine the difference between the (reduced) actions of the partner 
orbits. Apart from the difference associated to the energy offset already seen in the 
diagonal approximation this has a further contribution accounting for the change of 
action due to the changed connections in the encounters. For an encounter of two parts 
the latter contribution can be written as [35J [36, 137] 

A S = ^2 s k u k (46) 

k 


where the sum runs over pairs of associated stable and unstable coordinates. 
Generalisations to encounters involving more orbit parts are given in [IT]. As in 
the phase of (42) the action difference is divided by h, systematic contributions of 
constructively interfering orbits will have action differences of this scale and hence very 
close encounters. 

As a further ingredient one needs to determine the probability that encounters 
with given separations between the stretches arise in a periodic orbit/solution. For the 
corresponding formula we refer to PH- We stress that the only requirement for its 
validity is the existence of spectral gap. This is used to derive an ‘ergodic’ probability 
for different parts of an orbit to come close, as well as an estimate for the duration of 
an encounter, both of which enter the formula. 

With all ingredients for the single-particle treatment remaining valid, the 
contributions of all ‘diagrams’ of orbits differing in encounters (such as those displayed 
in Fig. [l| remain unchanged. These contributions involve powers of 1/e with the power 
increasing for more complex diagrams. The sum over all diagrams can be performed 
using the combinatorial techniques discussed in m- For systems without time-reversal 
invariance the contributions cancel leaving only the diagonal approximation. For time- 
reversal invariant systems one obtains a series expansion to all orders in 1/e. In either 
case the results agree with the non-oscillatory terms with coefficients c n in the random 
matrix prediction 

T 


R(e) ~ 1 + Re ^ ^ (c n + d n 


„2ze\ 


n =2 


(47) 


Here systems without time-reversal invariance (as described by the Gaussian Unitary 
Ensemble) have C 2 = —^ and = \ and all other coefficients vanish. For time-reversal 
invariant systems (as described by the Gaussian Orthogonal Ensemble) we have C 2 = —1, 
c n = (n ~ 3 2 -l n ~ 1) for n > 3, d 2 = d 3 = 0 and d n = for n > 4. We note that the 

treatment summarised here assumes that our system (in the reduced phase space) has 
no further symmetries as these would lead to additional pairs of orbits with similar or 
identical actions, such as mutually reflected orbits in a system with reflection symmetry. 
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The oscillatory terms and systems with discrete geometric symmetries will be discussed 
later. 

As in case of normal ordering the trace formula is modified to include the Solari- 
Kochetov (SK) phase, we have to check that this does not affect the results arising in 
the present approach. To do so, we note that the contributions to the double sum only 
involve differences between phases associated to closeby partner orbits. Due to -0 oc ^ 
the SK phase itself is of an order independent of h. As the SK phase is time-reversal 
invariant the phases of two partner orbits can only differ due to the slight changes inside 
encounters. However for the encounters relevant for spectral statistics the encounters 
become very close in the semiclassical limit, meaning that the difference between the 
SK phases becomes negligible. 


3 . 4 . Oscillatory contributions 


The random-matrix prediction (47) for R(e) also involves oscillatory contributions 

■ ,2 * e To access these contributions a more careful semiclassical 


proportional to Re^Ay^e 2 
approximation is needed 
from spectral determinants via 


In this approximation the level density is accessed 


d (B) = V fl det(E - H) 


vr drj det(i? + r] — H) 77=0 
An approximation for the spectral determinant on the level of the trace formula is 


(48) 


det(£ - H ) cx exp ( - / dE'tr(E' - H) 


\-l 


oc e 


-inN(E) 


5>r(-l) 


nr e iS r (E)/h . 


(49) 


Here the sum is taken over sets of orbits T with nr elements and cumulative reduced 
action Sf; the amplitude Fr depends on the stability and the Maslov indices of the 
contributing orbits and N(E) is the smooth approximation for the number of energy 
levels below E. However using the spectral determinant allows to incorporate further 
quantum mechanical information, in particular the fact that due to H being Hermitian 
det(A — H) has to be real for real arguments E. As shown in [38] this leads to an 
approximation for the spectral determinant where the contributions from sets of orbits 
with cumulative periods larger than half of the Heisenberg tine Th = 2nhd are replaced 
by the complex conjugate of the contributions from sets with cumulative periods below 
this threshold. This ‘Riemann-Siegel lookalike formula’ readily generalises to the many- 
body systems under consideration as its key ingredient, the semiclassical approximation 
for the trace of the resolvent (E — H) -1 , can be accessed from the trace of the propagator 


as above; essentially one only has to omit the restriction to the imaginary part in (13). 
Again periodicity is required only in the reduced phase space. 

Using the improved approximation of the level density as well as the same general 
ideas as outlined above it is possible to resolve oscillatory contributions as well. As each 
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level density brings in two determinants, evaluating the two-point correlation function 
then requires to study interference between quadruplets of (possibly empty) sets of 
orbits. Hence one also needs to take into account contributions where, say, after changing 
connections inside encounters an orbit is broken into two orbits with similar cumulative 
action. A treatment of these more involved correlations shows that for chaotic quantum 
systems R(e) fully agrees with the predictions from RMT jT8lj . 


3.5. Systems without a spectral gap 


Interestingly, there are chaotic systems that do not have a spectral gap but still satisfy 
some of the ingredients of our calculation, such as hyperbolicity which is required for the 
existence of orbit pairs differing in encounters. For these systems many of the techniques 
sketched here are applicable, but the final result of Wigner-Dyson statistics does not 
carry over as e.g. the sum rule (44) becomes invalid. An interesting question for future 


work is whether chaotic many-body systems without a spectral gap could be faithful 
to RMT ensembles that incorporate more information about the problem at hand and 
reduce to Wigner-Dyson statistics in important regimes, e.g. the embedded many-body 
ensembles mm or ensembles sensitive to the spacial structure of the problem. An 
example for a situation where the spacial structure is important is the continuum limit 
with diverging number of sites; additional orbit correlations relevant in this case were 
identified in 


4. Symmetries 

In a system with additional discrete symmetries one needs to consider the spectral 
statistics inside subspectra determined by the symmetry group. A prime example is 
the discrete (disorder-free) Bose-Hubbard model with periodic boundary conditions as 
introduced above. Its symmetry group is the dihedral group, consisting of the discrete 
translation ijjk —» t'Wi and its iterates; the reflection tfk —>• ipL-i-k] and combinations 
of translations and reflection that can be viewed here as reflection about a different 
centre. More formally the symmetry group for this model is the dihedral group Dl- 
Here D n stands for the dihedral group of order 2 n. The spectrum of the Bose-Hubbard 
model decomposes into subspectra labelled by the eigenvalues e 27r * K / L of the discrete 
translation operator where k — 0,..., L — 1 denotes the quasi-momentum. The spectra 
with k = 0 and (if integer) k = L /2 decompose further into components even and odd 
under reflection, whereas the remaining subspectra come in energy-degenerate pairs k, 
L — k related by reflection. Altogether we obtain L + 1 subspectra for odd L and L + 2 
subspectra for even L. A representation-theoretic justification of this decomposition will 
be given below and we refer to [42J for ways to implement the decomposition numerically. 

The level density associated to each subspectrum is obtained using a trace formula 
where the classical orbits, or in the present context the solutions of the Gross-Pitaevski 
equation, are required to be periodic (in the sense above) inside a fundamental domain 
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of the system [33]. (See also [44j.) For the case at hand this domain can be defined e.g. 
by converting if' k to ifk for a fixed choice of ip' 0 and then imposing certain conditions on 
if k . Demanding that Keifk is smallest for k = 0 guarantees that applying the translation 
operator to a point inside the fundamental domain leads to a point outside. The 
same applies to most other elements of the symmetry group but in order to guarantee 
that reflection about the 0-th site leads outside the fundamental domain we need an 
additional requirement such as Re^/h < ReipL-i- 

Each subspectrum a (not distinguishing between degenerate subspectra) can 
now be associated to an irreducible representation of the symmetry group, and the 
corresponding trace formula [43] has a form similar to (JTj) , 

d a (E) « d a (E) + ^Re £ X a (g P )A p e iS ^ h . (50) 

p 


The only difference from ([Tj) apart from the restriction to the fundamental domain is 
the additional factor Xa(g P )- Here g p is the group element relating the initial and final 
point of the orbit p if the orbit is considered in the full phase space as opposed to the 
fundamental domain. The character Xa(g P ) is the trace of the matrix representing g p in 
the representation a. As in [43] the derivation of (50) requires that a projection operator 


on the part of the Hilbert space associated to our subspectrum is inserted inside the 
trace taken over the Hilbert space; doing this in our many-particle calculations starting 
from (13) leads to exactly the same modifications as observed in [43] for single-particle 
systems. 


To apply (50) to the Bose-Hubbard model we need the irreducible representations 
of the dihedral group m- The representation of the translation operator must have an 
L-fold power I; its two-dimensional representations are therefore be given by 



(51) 


The reflection operator is then represented by 



(52) 


and the representations of all other group elements can be found using that the matrix 
representation of a product of symmetry operations must be the product of the matrix 
representations. The general theory of symmetries in quantum mechanics (see e.g. [46]) 
implies that the energy eigenstates associated to this representation can be grouped 
into pairs (representable as vectors) with identical energy. One can show that for the 
two-dimensional representation at hand this leads precisely to the energy-degenerate 
subspectra associated to eigenvalues of the translation operator e 2mK ^ L (k 0, k) as 
mentioned above. 

However if we have k = 0 or, assuming even L, k = ^ the matrix (51) becomes 
diagonal and the two-dimensional representations defined above become reducible. 
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Instead of using these representations the associated subspectra therefore have to be 
described by one-dimensional representations. These represent the translation operator 
either by 1 or in case of even L alternatively by -1, whereas the reflection operator can be 
represented by either 1 or —1 regardless of L. Again these spectra precisely correspond 
to what we said above about the cases k — 0, 

In general the consideration of discrete symmetries can change the appropriate 
RMT ensemble compared to the one expected based on the time-reversal properties 
alone 03 • However using Eq. fl50| ) one can show that this does not happen for 

As the aforementioned 


subspectra associated to representations by real matrices 
representations of the dihedral group are real the spectral statistics of all its subspectra 
is thus in line with the GOE as observed numerically in as well as in the next 
section. 


5. Numerical results 


To support our results numerically and complement the previous numerical studies we 
have performed a numerical analysis of the chaotic properties of both the quantum and 
the classical Bose-Hubbard model. Our results support Wigner-Dyson statistics as well 
as mostly classical chaotic dynamics in the regime where the hopping and interaction 
terms are comparable. 

For the quantum model we are interested in the spectral statistics as discussed 
above. We use a statistical observable slightly more convenient for computations 
than R(e), namely the normalised distribution P{r ) of ratios between subsequent level 
spacings; if the ordered quantum levels are denoted by E n , these ratios are given by 

r n = " +1 —-. The ratio distribution is especially suited for our purposes as there 

En E n —i 

is no requirement to explicitly evaluate the average level density. A random matrix 
prediction for P(r) was obtained in [50J by considering 3x3 random matrices; for the 
GOE it reads 


Prmt{x) 


27 r + r 2 
8 (1 + r + r 2 ) 5 / 2 


(53) 


In particular for r —>■ 0 one can see that P(r) oc r, which is due to level repulsion. 
Notice that for large r the distribution has a fat tail in contrast to the level spacing 
distribution: P(r) oc r~ 3 for r 1. As for the density of nearest-neighbour spacings the 
results for large matrices are expected to be very similar in value but analytically more 
complicated. For comparison, we mention that the ratio distribution for an integrable 
system (assuming that the energy levels are independent) is given by 


P, 


Poisson 


(r) = 


(1 + rf 


(54) 


We determined the spectrum of the quantum Bose-Hubbard model via exact 
diagonalisation. For each of the subspectra described in section [| we computed the 
histogram of r to get a numerical estimate of the ratio distribution P num (r). Examples 
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Figure 2. Numerical estimate of the ratio distribution between neighbouring levels 
for the quantum Bose Hubbard model with L = 5 and N = 25. The dashed black line 
stands for the RMT prediction (53). The dotted-dashed red line stands for the Poisson 
distribution ( |54| ). Green circles: UN/J = 0.125. Blue squares: UN/J = 5. Orange 
triangles: UN/J = 250. 


of such numerical histograms are displayed in Fig. [2j Then we determined the L 1 norm 
(i.e. the integral over the absolute value) of the difference P num (r ) — Prmt(j')- Finally 
this difference was averaged over all different subspectra. Fig. [3] top shows the norm 
for the case L — 5 and N = 15 as well as iV = 25. We obtain good agreement for 
the case that the interaction term (of order UN 2 ) is comparable to or slightly larger 
than the hopping term (of order JN). The agreement improves as N increases, which 
is reassuring as our theoretical derivation of Wigner-Dyson statistics holds in the limit 
N — y oo (or equivalently in the semiclassical limit). 

For the classical Bose-Hubbard model we introduce a numerically determined 
measure of ergodicity, which indicates whether the classical limit is mostly chaotic. 
To determine this measure we considered the Hamiltonian 




j 

2 


L—l jj L—l 

+ Vk^k+l) + j 


where we identified ipi = i/)q and f>* L = ip/). The classical dynamics is given by 


ihxj> = 


dH 

dt/j* 


ihif 


dH 
dip ’ 


(55) 


(56) 
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UN/J 


Figure 3. Top: Comparison between the ratio distribution of neighbouring levels 
between the Bose-Hubbard model and RMT for N = 25 (black circles) and TV = 15 
(red squares), and L = 5. The difference is estimated via the L 1 norm of the deviation 
from RMT as a function of UN/J. Bottom: Numerical estimate of the degree of 
ergodicity of the classical dynamics of Bose-Hubbard model with N = 25 and L = 5 as 
a function of UN/J. 1 indicates full ergodicity. The precise definition of our estimate 
is given in the text. 


If an initial point {if q, i/’q) is given in phase space, Eq. (56) allows to compute the unique 
trajectory (i/>*(t), ip(t)) such that (i/>*(0), i/>(0)) = (i/>o, ^o)- It may be useful to add 
that as each point is part of a classical trajectory the point (i/>g, ?/> 0 ) is parameterised by 
2 L real numbers as the coefficients of i/>q are the complex conjugates of the coefficients 
of i/v 

As the energy E and the particle number N are conserved, the trajectory can only 
access a subspace in phase space, denoted by Ai n,e, which is the constant energy surface 
associated to E in the Fock space with N particles. Numerically a triangulation of Ai n,e 
is performed by picking random points on this surface. These points are generated using 
Halton sequences, which are commonly used to explore high dimensional spaces. First a 
random point is sampled on the hypersphere | |?/>| | J = N. If its energy coincides with E 
it is kept and stored as a point in Ai n,e■ Otherwise we continue the sampling. In order 
to cover substantially M.n,e for typical values of E, UN/J we chose to pick 2 x 10 6 
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random points. Then the mean distance between any two of these points is calculated, 
let us call it 5. This gives a typical distance scale for the constant energy surface at a 
given energy E. This typical distance is used to define balls of radius 8 /10 around each 
points. After having determined a cover of M.n,e the trajectory starting from (if) q, t/? 0 ) 
is computed for a sequence of increasing (final) times. We took large enough final values 
of the propagation time so that the number of visited balls does not vary significantly 
after this time. 

Finally a measure for the ergodicity of a trajectory is defined by the ratio between 
the number of visited balls and the total number of balls. In order to have generic 
values, this ergodicity measure is averaged over different choices of initial conditions 
(i/>o, t/> 0 ). Our results are displayed in Fig. [3] bottom. By definition our measure is a 
real number between 0 and 1. 1 indicates perfect ergodicity. The closer it is to 1 the 
more ergodic the trajectory is. While not a rigorous check of the conditions for random 
matrix statistics, our measure indicates whether the system is mostly ergodic or has a 
substantially mixed phase space with large stability islands. Fig. [3] bottom shows that, 
for the same range of UN/J , the classical dynamics is predominantly ergodic and the 
quantum system agrees closely with RMT prediction. 

Our results are compatible with those obtained through a different approach and a 
modified version of the model in [ST] . There the authors chose a different normalisation 
of t/j and claimed that their version of the classical Bose-Hubbard model becomes more 
and more classical when increasing U/ J. This is line with the left part of Fig. [3] bottom. 
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Appendix A. Perturbation of the stability matrix and the Hessian 


In this appendix we want to explicitly show how perturbation theory can be used to 
deal with the vanishing eigenvalues of the stability matrix and the Hessian considered in 
subsection 2.5, associated to k — 0 and k — 1. To avoid some of the complexity of [25] 
we specifically consider small perturbations 8H that depend only on 8' k , which does not 
appear as a parameter of the unperturbed Hamiltonian. In order to avoid ambiguities 
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we choose SH(9' k ) periodic in 9' k with the period coinciding with the range of values 
covered by 9' k . 

We first investigate how the perturbation modifies the vanishing eigenvalue of the 
Hessian %. The corresponding eigenvector is associated to a simultaneous shift of 9'- k 
for our k at all time steps j. It is thus convenient to rewrite the Hessian in terms 
of derivatives w.r.t. Ij k , 6'- k ; this also absorbs the divisor li in (21) and it has the 


advantage that the Hessian becomes Hermitian. The perturbation of the Hessian d'H 

. . , ... . .... d 2 6H(6') 

has as its only non-vanishing matrix elements the derivatives- — 




-t. In the same 


coordinate system the eigenvector e associated to the vanishing eigenvalue has identical 
entries for the J components associated to the 9[ jk associated to our k and the remaining 
components are zero. First-order perturbation theory then yields a perturbed eigenvalue 


e ■ ddie 


e ■ e 


-~y 

j z - 


d^H(9' k ) 

d0% 


(A.l) 


r(U 


taken as e in subsection | 

To study the effect on the stability matrix we first consider the matrix ATJ ' mapping 
a deviation of (#' fc , Ij k ) to the resulting deviation after a time step of size r. We obtain, 


e.g. by converting Eq. (27) to the present system of coordinates, 


= 


i , sn h0 2 j j ) t 

\ ar/ 1 


d 2 H(0j ,Ij 


't'lr.o .1 ■ 
as % 


■T 


T 


1 - 


ai % 

d 2 H(ej,I j ) 

a Vik 


T 


(A.2) 


Here the lower left entry vanishes for the unperturbed Hamiltonian but is moved away 
from zero by the perturbation SH(9 k ). As a consequence of this perturbation the product 


given in (32) 


M {k) = 


dA9' k 

dl’ 


(A.3) 


then receives the lower left entry 


-E 


80% 


t = ej 


(A.4) 


as desired. In contrast the changes of the other entries do not affect the determinant to 
linear order in the perturbation. 


Appendix B. Weyl ordering 

Finally we want to show explicitly that in line with [271 [231. 25l [9J [33] the SK phase 
does not arise in case of Weyl ordering. Using the Hamiltonian in Weyl ordering the 
short-time propagator can be written as [33| 














Spectral statistics of chaotic many-body systems 


26 


~ 2i J d \ w P w *j\ ex P ( - 2| Wj\ 2 + 2■ Wj + 2^ j _ l ■ w* 

-l^il 2 / 2 - l^-il7 2 - <7 ' *l>j -1 - ^H(w),Wj)T} (B.l) 

involving an integration over a pair of complex conjugate variables Wj, w* that were not 
required in normal ordering. The action in the discretised coherent-state path integral 
then turns into [33] 

R = (“( _ 2 K '| 2 + 2 ^j • W 3 + 2 7 j-l • W *j - l^jl 2 / 2 - l^j-ll 2 / 2 


(B.2) 


^r^j- 1 ) - 


and the path integral runs over w v w* as well as if if*. The action becomes stationary 
under variation of -*/?•, ?/>* if 

7, + Vb-1 


Wj = 


(B.3) 


and its complex conjugate hold. This suggests that in the limit of fine discretisation 
the Wj should be interpreted as macroscopic wave functions halfway between two 
discretisation steps. Furthermore stationarity w.r.t. variation of Wj,w* implies 


ih- 


Wj 


*l>i -i dH(w*,Wj ) 


r/2 


dw* 


(B.4) 


and its complex conjugate, which is the appropriate discretisation of the nonlinear 
Schrodinger equation ihxj) = for the time interval r/2. After eliminating Wj with 
the help of (B.3) a short calculation shows that the combination of macroscopic wave 


functions in the action (B.2) agrees with the simpler form in (12) 


We now have to evaluate the Hessian of j and its determinant. In analogy to 


subsection 2.4 we start with the Hessian containing all second derivatives. If we first 
write the derivatives w.r.t. all w*,Wj and then those w.r.t. all the Hessian 

_ 2E 

T T ] where D is block-diagonal with the blocks 

-2 E E + E 1 


assumes the form i 


Dj — 2I2 l + 


it d 2 H(w*,Wj) 
h d(w*,Wj ) 2 

and the orthogonal matrix E is given by 

/0 I L \ 

E = 


(B.5) 




II 


(B.6) 
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det i 


D -2 E 

-2 E t E + E t 

= det(2 E t ) det(2 E - D(2E T )~\-E - E T )) 
= det (D + DE 2 - 4 E) 

( D 0 + AF D 0 + 4 F t 

D\ +4 F D\ + 4 F t 


= det 


VA7-1+4F 2 


\ 

Dj-\ + AF J 


with F = 


0 If 
0 0 


. Here the factor i dropped out because the dimension of the matrix 


is a multiple of 4. As the final matrix is of block tridiagonal form it can be simplified 
with the help of (24). Using 


det(.Dj + AF r ) =4 L + 0(r 2 ) 

(Dj + ^F T )-\D j + AF) = Mj + 0 (t 2 ) 


(where Mj is defined as in (27) but with Wj , w* taking the role of xfj , ip *) and again 
ignoring terms of quadratic and higher order in r in the entries we obtain 


det 



2 E \ 
-E - E t ) 


A jl det(M-I 2L ). 


(B.7) 


As usual the stationary-phase approximation requires the inverse square root of this 
determinant, containing a factor 2~ JL . This factor is cancelled by the 2^ L in the 
integration measure for Wj , w*. As anticipated, we thus obtain the same result as with 
normal ordering apart from the SK phase. The remaining steps, in particular taking 
into account the conservation laws, carry over directly from the case of normal ordering 
discussed in the main text. 
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